library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
library(ggpubr)
library(dbscan)
##
## Attaching package: 'dbscan'
## The following object is masked from 'package:stats':
##
## as.dendrogram
library(corrplot)
## corrplot 0.94 loaded
library(corrr)
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
theme_set(theme_minimal())
Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.
df_full <- readRDS("~/BioStat/Прод_виз_данных/very_low_birthweight.RDS")
df_cleaned1 <- df_full %>%
select(where(~ sum(is.na(.)) <= 100)) %>%
filter(if_all(everything(), ~ !is.na(.)))
summary(df_full)
## birth exit hospstay lowph
## Min. :81.51 Min. :68.53 Min. :-6574.00 Min. :6.530
## 1st Qu.:83.52 1st Qu.:83.58 1st Qu.: 16.00 1st Qu.:7.130
## Median :84.90 Median :84.96 Median : 37.00 Median :7.210
## Mean :84.75 Mean :84.84 Mean : 40.36 Mean :7.202
## 3rd Qu.:86.07 3rd Qu.:86.17 3rd Qu.: 62.00 3rd Qu.:7.310
## Max. :87.48 Max. :96.87 Max. : 3668.00 Max. :7.550
## NA's :21 NA's :31 NA's :31 NA's :62
## pltct race bwt gest
## Min. : 16.0 black :369 Min. : 400 Min. :22.00
## 1st Qu.:143.0 native American: 16 1st Qu.: 900 1st Qu.:27.00
## Median :202.0 oriental : 4 Median :1120 Median :29.00
## Mean :201.6 white :257 Mean :1094 Mean :28.87
## 3rd Qu.:252.0 NA's : 25 3rd Qu.:1310 3rd Qu.:31.00
## Max. :571.0 Max. :1580 Max. :40.00
## NA's :70 NA's :2 NA's :4
## inout twn lol magsulf
## born at Duke:547 Min. :0.0000 Min. : 0.000 Min. :0.0000
## transported :121 1st Qu.:0.0000 1st Qu.: 0.000 1st Qu.:0.0000
## NA's : 3 Median :0.0000 Median : 3.500 Median :0.0000
## Mean :0.2074 Mean : 8.438 Mean :0.1344
## 3rd Qu.:0.0000 3rd Qu.: 9.000 3rd Qu.:0.0000
## Max. :1.0000 Max. :192.000 Max. :1.0000
## NA's :20 NA's :381 NA's :247
## meth toc delivery apg1
## Min. :0.0000 Min. :0.0000 abdominal:314 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.0000 vaginal :335 1st Qu.:2.000
## Median :0.0000 Median :0.0000 NA's : 22 Median :5.000
## Mean :0.4372 Mean :0.2248 Mean :4.903
## 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:7.000
## Max. :1.0000 Max. :1.0000 Max. :9.000
## NA's :106 NA's :106 NA's :34
## vent pneumo pda cld
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :1.0000 Median :0.0000 Median :0.0000 Median :0.0000
## Mean :0.5803 Mean :0.1969 Mean :0.2087 Mean :0.2694
## 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
## NA's :30 NA's :26 NA's :29 NA's :66
## pvh ivh ipe year sex
## absent :360 absent :442 absent :472 Min. :81.51 female:320
## definite:125 definite: 75 definite: 38 1st Qu.:83.52 male :330
## possible: 41 possible: 10 possible: 17 Median :84.91 NA's : 21
## NA's :145 NA's :144 NA's :144 Mean :84.76
## 3rd Qu.:86.07
## Max. :87.48
## NA's :21
## dead
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.2146
## 3rd Qu.:0.0000
## Max. :1.0000
##
str(df_full)
## 'data.frame': 671 obs. of 26 variables:
## $ birth : num 81.5 81.5 81.6 81.6 81.6 ...
## $ exit : num 81.6 81.5 81.6 81.7 81.6 ...
## $ hospstay: int 34 9 -2 40 2 62 32 NA NA 28 ...
## $ lowph : num NA 7.25 7.06 7.25 6.97 ...
## $ pltct : int 100 244 114 182 54 NA 282 NA NA 153 ...
## $ race : Factor w/ 4 levels "black","native American",..: 4 4 1 1 1 4 1 NA NA 1 ...
## $ bwt : int 1250 1370 620 1480 925 940 1255 600 700 1350 ...
## $ gest : num 35 32 23 32 28 28 29.5 26 24 34 ...
## $ inout : Factor w/ 2 levels "born at Duke",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ twn : int 0 0 0 0 0 0 0 NA NA 0 ...
## $ lol : int NA NA NA NA NA NA NA NA NA NA ...
## $ magsulf : int NA NA NA NA NA NA NA NA NA NA ...
## $ meth : int 0 1 0 1 0 1 1 NA NA 1 ...
## $ toc : int 0 0 1 0 0 0 0 NA NA 0 ...
## $ delivery: Factor w/ 2 levels "abdominal","vaginal": 1 1 2 2 1 1 2 NA NA 1 ...
## $ apg1 : int 8 7 1 8 5 8 9 NA NA 4 ...
## $ vent : int 0 0 1 0 1 1 0 NA NA 0 ...
## $ pneumo : int 0 0 0 0 1 0 0 NA NA 0 ...
## $ pda : int 0 0 0 0 0 0 0 NA NA 0 ...
## $ cld : int 0 0 NA 0 0 0 0 NA NA 0 ...
## $ pvh : Factor w/ 3 levels "absent","definite",..: NA NA NA NA 2 1 NA NA 1 NA ...
## $ ivh : Factor w/ 3 levels "absent","definite",..: NA NA NA NA 2 1 NA NA 1 NA ...
## $ ipe : Factor w/ 3 levels "absent","definite",..: NA NA NA NA NA 1 NA NA 1 NA ...
## $ year : num 81.5 81.5 81.6 81.6 81.6 ...
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 2 1 1 1 NA NA 1 ...
## $ dead : int 0 0 1 0 1 0 0 1 1 0 ...
Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.
df_cleaned2 <- df_cleaned1 %>%
mutate(
twn = as.factor(twn),
vent = as.factor(vent),
pneumo = as.factor(pneumo),
pda = as.factor(pda),
cld = as.factor(cld),
dead = as.factor(dead),
apg1 = as.factor(apg1)
) %>%
mutate(across(where(is.numeric), ~ ifelse(
. < quantile(., 0.25, na.rm = TRUE) - 1.5 * IQR(., na.rm = TRUE) |
. > quantile(., 0.75, na.rm = TRUE) + 1.5 * IQR(., na.rm = TRUE),
NA, .))) %>%
filter(if_all(everything(), ~ !is.na(.))) %>%
mutate(ID = row_number()) %>%
relocate(ID, .before = 1) %>%
filter(hospstay >= 0)
df2 <- df_cleaned2 %>%
select(where(is.numeric), -ID) %>%
pivot_longer(cols = everything(), names_to = "variable", values_to = "value")
ggplot(df2, aes(x = value)) +
geom_density(fill = "lightblue", alpha = 0.5) +
facet_wrap(~ variable, scales = "free") +
labs(title = "Density plots for numeric data", x = "", y = "Density")
plot1 <- ggplot(df_cleaned2, aes(x = pltct, fill = factor(inout))) +
geom_density(alpha = 0.5) +
labs(title = "Distribution density for lowph by variable 'inout'",
x = "lowph", y = "density", fill = "inout")
plot2 <- ggplot(df_cleaned2, aes(x = bwt, fill = factor(inout))) +
geom_density(alpha = 0.5) +
labs(title = "Distribution density for btw by variable 'inout'",
x = "bwt", y = "density", fill = "inout")
combined_plot <- ggarrange(plot1, plot2,
ncol = 2, # Количество столбцов
nrow = 1, # Количество строк
common.legend = TRUE, # Общая легенда для графиков
legend = "bottom" # Размещение легенды внизу
)
combined_plot
df_long <- df_cleaned2 %>%
select(ID, where(is.numeric)) %>%
pivot_longer(cols = -ID, names_to = "variable", values_to = "value")
ggplot(df_long, aes(x = ID, y = value)) +
geom_point(fill = "lightblue", alpha = 0.5) +
facet_wrap(~ variable, scales = "free")
Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?
wilcox_test_result <- df_cleaned2 %>%
wilcox_test(lowph ~ inout) %>%
mutate(y.position = 7.6) %>%
print()
## # A tibble: 1 × 8
## .y. group1 group2 n1 n2 statistic p y.position
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 lowph born at Duke transported 412 71 19947 0.000000957 7.6
ggboxplot(df_cleaned2, x = "inout", y = "lowph",
color = "inout", palette = "jco") +
stat_pvalue_manual(wilcox_test_result, label = "p", tip.length = 0.01) +
labs(title = "Wilcoxon Rank Sum Test", y = "lowph", x = "inout")
P-value из Wilcoxon теста (5.77e-08) значительно меньше 0.05, следовательно есть статистически значимая разница в значениях lowph между двумя группами: born at Duke и transported. Если более низкие значения lowph ассоциированы с худшей выживаемостью, можно предположить, что младенцы из born at Duke имеют лучшие шансы на выживаемость по сравнению с группой transported.
Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций. Постройте иерархическую кластеризацию на этом датафрейме.
df4 <- df_cleaned2 %>%
select(
where(is.numeric),
-c(ID, birth, year, exit)
)
cor_matrix <- cor(df4, use = "complete.obs")
print(cor_matrix)
## hospstay lowph pltct bwt gest
## hospstay 1.0000000 -0.2628111 -0.13751643 -0.5222201 -0.43391064
## lowph -0.2628111 1.0000000 0.29279438 0.3103002 0.35928694
## pltct -0.1375164 0.2927944 1.00000000 0.2180920 0.02837645
## bwt -0.5222201 0.3103002 0.21809200 1.0000000 0.67667411
## gest -0.4339106 0.3592869 0.02837645 0.6766741 1.00000000
plot4.1 <- df4 %>%
cor() %>%
corrplot(
order = 'hclust'
)
plot4.2 <- df4 %>%
cor() %>%
network_plot(min_cor = .0)
print(plot4.1)
## $corr
## hospstay pltct lowph bwt gest
## hospstay 1.0000000 -0.13751643 -0.2628111 -0.5222201 -0.43391064
## pltct -0.1375164 1.00000000 0.2927944 0.2180920 0.02837645
## lowph -0.2628111 0.29279438 1.0000000 0.3103002 0.35928694
## bwt -0.5222201 0.21809200 0.3103002 1.0000000 0.67667411
## gest -0.4339106 0.02837645 0.3592869 0.6766741 1.00000000
##
## $corrPos
## xName yName x y corr
## 1 hospstay hospstay 1 5 1.00000000
## 2 hospstay pltct 1 4 -0.13751643
## 3 hospstay lowph 1 3 -0.26281113
## 4 hospstay bwt 1 2 -0.52222009
## 5 hospstay gest 1 1 -0.43391064
## 6 pltct hospstay 2 5 -0.13751643
## 7 pltct pltct 2 4 1.00000000
## 8 pltct lowph 2 3 0.29279438
## 9 pltct bwt 2 2 0.21809200
## 10 pltct gest 2 1 0.02837645
## 11 lowph hospstay 3 5 -0.26281113
## 12 lowph pltct 3 4 0.29279438
## 13 lowph lowph 3 3 1.00000000
## 14 lowph bwt 3 2 0.31030019
## 15 lowph gest 3 1 0.35928694
## 16 bwt hospstay 4 5 -0.52222009
## 17 bwt pltct 4 4 0.21809200
## 18 bwt lowph 4 3 0.31030019
## 19 bwt bwt 4 2 1.00000000
## 20 bwt gest 4 1 0.67667411
## 21 gest hospstay 5 5 -0.43391064
## 22 gest pltct 5 4 0.02837645
## 23 gest lowph 5 3 0.35928694
## 24 gest bwt 5 2 0.67667411
## 25 gest gest 5 1 1.00000000
##
## $arg
## $arg$type
## [1] "full"
print(plot4.2)
Постройте иерархическую кластеризацию на этом датафрейме.
df_scaled <- scale(df4)
rownames(df_scaled) <- df_cleaned2$ID
dist_matrix <- dist(df_scaled, method = "euclidean")
hc <- hclust(dist_matrix, method = "ward.D2")
fviz_dend(hc,
cex = 0.6)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#fviz_dend(hc,
#k = 3,
#k_colors = "jco",
#type = "phylogenic",
#repel = TRUE) # Избежать наслоения лейблов
Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.
library(pheatmap)
pheatmap(df_scaled,
clustering_method = "ward.D2", # Метод кластеризации
scale = "none", # Данные уже нормализованы
color = colorRampPalette(c("navy", "white", "firebrick"))(50), # Цветовая палитра
main = "Heatmap with Hierarchical Clustering")
Субъекты с преимущественно низкими значениями hospstay могут представлять наблюдения с большими значениями bwt, gest и lowph и наоборот.
Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?
df.pca <- prcomp(df4,
scale = T)
summary(df.pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.5461 1.0294 0.8393 0.7541 0.5263
## Proportion of Variance 0.4781 0.2119 0.1409 0.1137 0.0554
## Cumulative Proportion 0.4781 0.6900 0.8309 0.9446 1.0000
fviz_eig(df.pca,
addlabels = T,
ylim = c(0, 50))
fviz_pca_var(df.pca, col.var = "contrib")
Первая компонента объясняет 47.8% дисперсии данных, первые две
компоненты объясняют 69% дисперсии данных, можно сделать вывод что их
можно использовать для дальнейшего анализа. Шкалирование применять
нужно, так как данные имеют разные единицы измерения и разную
дисперсию.
Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.
library(ggbiplot)
biplot <- ggbiplot(df.pca,
scale=0, alpha = 0.7,
groups = as.factor(df_cleaned2$dead))+
scale_color_discrete(
name = "Status",
labels = c("0" = "Alive", "1" = "Dead")
) +
labs(title = "PCA Biplot colored by 'dead'")
biplot
Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
pca_data <- as.data.frame(df.pca$x) # Извлекаем главные компоненты
pca_data$ID <- rownames(df_cleaned2) # Добавляем ID пациентов
pca_data$dead <- as.factor(df_cleaned2$dead)
pca_data$dead <- factor(pca_data$dead, levels = c("0", "1"), labels = c("Alive", "Dead"))
plot9 <- ggplot(pca_data, aes(x = PC1, y = PC2, color = dead, text = ID)) +
geom_point(size = 2, alpha = 0.7) +
labs(title = "PCA Biplot colored by 'dead'",
x = "PC1",
y = "PC2",
)+
scale_color_discrete(
name = "Status",
labels = c("0" = "Alive", "1" = "Dead"))
# Преобразование в интерактивный график с помощью plotly
plotly_plot <- ggplotly(plot9, tooltip = "text") # 'text' отображает ID пациента
plotly_plot
Дайте содержательную интерпретацию PCA анализу. Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?
Вообще (PCA) используется для уменьшения размерности данных, выявления главных направлений вариации в данных и для упрощения визуализации данных. Как было сказано ранее, первая компонента объясняет 47.8% дисперсии данных, первые две компоненты объясняют 69% дисперсии данных, что является хорошим показателем для дальнейшего анализа. Самые длинные стрелки у переменных hospstay и gest, они вкладывают больше всего информации в объяснении вариации, также видно что они направлены в одну сторону, что свидетельствует об их корреляции. Вклад hosptstay, bwt и gest сильнее в компоненту PC1, pltct и lowph влияют на PC2. Из графика видно, что точки с Alive и Dead значениями не образуют кластеров, что указывает отсутствие тенденции разделения выживших и умерших на основе первых двух главных компонент.
PCA основывается исключительно на входных переменных (hosptstay, bwt, gest, lowph, pltct) и их вариации. Переменная dead используется только для визуализации и никак не влияет на формирование главных компонент. Чтобы проверить, связаны ли переменные с выживаемостью корректнее использовать модели, такие как логистическая регрессия.